Nonlinear Matter Wave Dynamics with a Chaotic Potential 
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We consider the case of a cubic nonlinear Schrodinger equation with an additional chaotic poten- 
tial, in the sense that such a potential produces chaotic dynamics in classical mechanics. We derive 
and describe an appropriate semiclassical limit to such a nonlinear Schrodinger equation, using a 
semiclassical interpretation of the Wigner function, and relate this to the hydrodynamic limit of 
the Gross-Pitaevskii equation used in the context of Bose-Einstein condensation. We investigate a 
specific example of a Gross-Pitaevskii equation with such a chaotic potential: the one-dimensional 
^\ • delta-kicked harmonic oscillator, and its semiclassical limit. We explore the feasibility of experimen- 

Q\ ' tal realization of such a system in a Bose-Einstein condensate experiment, giving a concrete proposal 

ON ■ of how to implement such a configuration, and considering the problem of condensate depletion. 
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O ! I. INTRODUCTION 

(N 

Chaos in classical Hamiltonian systems, most simply thought of as the extreme sensitivity of trajectories in phase 
| space to initial conditions, making long term predictions extremely difficult, is by now broadly understood jj],^]]. 
More recently, the field of quantum chaos, for our purpose meaning the study of quantum mechanical equivalents of 
classical chaotic systems, has been the subject of much investigation w)$i. From this it does seem that the dynamics 
of quantum mechanical systems can be divided into regular and irregular subsets, with distinct differences between the 
two, just as is the case in classical mechanics. For example, due to the unitarity of the evolution of the state vector, 
there can be no equivalent of sensitivity to initial conditions in the Hilbcrt space, but there appears to be an equivalent 
sensitivity to perturbation which distinguishes quantum chaotic motion Q . A certain amount of understanding has 
thus been achieved, although there are still unresolved problems, in particular how to extract classical chaos from 
quantum mechanics j^j. 

Quantum dynamics are determined by the Schrodinger equation. A seemingly natural extension is to ask what 
happens when we take a quantum chaotic Schrodinger equation, and add some kind of nonlinearity. This is something 
which has been much less studied and is certainly of more than academic interest; such equations do appear in 
nature, for example the Gross-Pitaevskii equation in the field of Bose-Einstein condensation , and also in the field 
of nonlinear optics ||] . There are thus experimentally accessible systems in which such chaotic effects may manifest 
themselves. 

It is also interesting to note that just as in the case of quantum mechanics, where if one takes the limit h — ► 
one expects to regain classical dynamics, one can also carry out this limit for nonlinear Schrodinger equations. This 
produces equations reminiscent of classical hydrodynamics, an interpretation also extensively used in the theoretical 
study of Bose-Einstein condensates jl0| . This interconnection of different kinds of dynamics is displayed schematically 
in Fig. ®. 

We will firstly be concerned with effective "single particle" systems. That is to say, where one takes a linear single 
particle Schrodinger equation with a potential which is known to produce chaotic dynamics in Hamilton's equations of 
motion, and adds a nonlinearity to it. The Gross-Pitaevskii equation (for example) describes the collective dynamics 
of huge numbers of particles, but may nevertheless be thought of as an effective single particle wave equation. We 
later consider corrections to this interpretation, taking into account more fully the many body dynamics. 

It is thus of general interest to determine how effects of classical chaos and quantum chaos manifest themselves 
in the dynamics of nonlinear Schrodinger equations, and to what extent the dynamics of the nonlinear Schrodinger 
equation can be explained by motion in the hydrodynamic limit, as determined by the hydrodynamic equations. 



II. GENERALITIES 
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A. Gross-Pitaevskii Equation 



In this paper we will consider explicitly only one dimensional systems, although the analytic results presented 
can easily be generalized to two or three spatial dimensions. To simplify things further, we consider only the cubic 
nonlinearity explicitly, the simplest nonlinearity possible, resulting in the one dimensional Gross-Pitaevskii equation, 
well known in the context of Bose-Einstein condensation: 

d ti 2 d 2 

where tp(x,t) is the wavefunction and u the strength of the nonlinearity. Again, the analytic results here can easily 
be generalized to more complicated nonlincarities. Such a simplified system demonstrates all the main features of a 
nonlinear Schrodinger equation, and is perfectly adequate for illustrative purposes. This kind of simplified system is 
in fact experimentally accessible, for example in a Bose-Einstein condensate experiment, as will be shown in Section 



B. Hydrodynamic Equations 



It is tempting to think of the hydrodynamic equations as the semiclassical limit of the Gross-Pitaevskii equation. 
This turns out to be not quite so, as will be shown in Sec. [II. We nevertheless sketch out the standard derivation of 



the hydrodynamic equations, in order to set notation, and so that later we can point out the differences between the 
hydrodynamic limit and the genuine semiclassical limit, which we will derive using Wigner functions. 

We rewrite the Gross-Pitaevskii equation Eq. (|l|) using the density p and a momentum field P , defined in terms 
of the wavefunction ip = ^fpe lS l h as: 
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The resulting equation of motion for the density is 

Before moving to the equation of motion for P, we first consider the equation for S, which is 
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The equation of motion for the momentum field P is exactly the spatial derivative of Eq. (|^): 
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Taking the hydrodynamic limit J8|,[l0| consists of abandoning the term in Eq. (||) proportional to ti 2 , generally 
justified by claiming that the density p is sufficiently smooth for its derivatives to be insignificant, resulting in 



dt dx 
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(7) 



Clearly, to get Eq. (|7|), we have discarded all quantum character of the Gross-Pitaevskii equation. Also note that 
if the corresponding term is abandoned in Eq. (^|) in the case where u — and V is time independent, we get the 
Hamilton- Jacobi equation for a single particle in the potential V, with the interpretation that dS/dx is the canonically 
conjugate momentum to the coordinate x [ p"2| . 

This seems to indicate that the hydrodynamic equations Eqs. might be an equivalent "classical" limit to the 
Gross-Pitaevskii equation with finite u [Eq. ([[])]. As previously stated however, this turns out to be not quite so, as 
we shall soon see. 
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III. WIGNER FUNCTION DYNAMICS 



A. Expansion in h 



We wish to carry out a consistent expansion of Eq. (|l|) around h, in order to clearly separate classical from 
quantum dynamics, and to provide order by order corrections. We will do this by considering the dynamics of the 
Wigncr function W , which is exactly equivalent to the wavefunction ip, in the sense that all information about the 
wavefunction is contained within its Wigner representation. 

We define the Wigner function (for a pure state) as 

W{x,p) = — / d Te ~ ipT/h (p*(x - r/2)(p(x + r/2). (8) 

It is well known that the dynamics of the Wigner function of a single particle to lowest order give simply the classical 
Liouville equation of a distribution of noninteracting particles || . The exact expression to all orders in H for the time 
evolution of the Wigner function W is given by: 

°-W = f J=£- ( h X ^H^W - —H—W (9) 
dt f^ o (2s + 1)\ \2 J dx 2 ^ 1 9p 2s +! dp dx K > 

where H is the single particle classical Hamiltonian function. How to obtain this expression is sketched in Appendix 
[A|. Setting h — we see we do indeed get the classical Liouville equation 

d d d d d 

dt^ dx^ dp^ dp^ dx^' 

so long as the initial Wigner function can in fact be interpreted as a classical probability density (i.e. is non-negative). 
If we have as a classical Liouville density a delta distribution, W(x,p) = S(x — xq)S(p — po), we regain classical point 
dynamics. One can think of a point particle being regained from quantum mechanics if we have a coherent state 
centred at x = xq and p — po and let h — > 0, causing the Wigner function to tend to just such a delta distribution. 

It is worth mentioning that although we talk blithely about letting h tend to zero, this is in fact physically 
meaningless. As H is a constant, we must in fact expand around some scaling parameter to do with the characteristic 
action scales of the problem at hand, such that at some point the quantum corrections should be completely dominated, 
at least for some characteristic time ||. Generally some appropriate parameter presents itself, as will be shown in 
the model we present in Section |V|, and expansions where it is stated that the limit /j, — ^ is explored should be 
interpreted in this manner. 

What we now wish to do is to take an equivalent limit to that presented in Eqs. (|^,|l0|) for the Gross-Pitaevskii 
equation, with the object of getting some kind of Liouville equation with the nonlinearity taken into account. The 
full expansion of the Wigner function dynamics governed by Eq. (]]]) in terms of h turns out to be: 

8 B 8 ^ (-'\] s fh\ 2s 8 2s+1 8 2s+1 

where we have the density 

p(x) = / dp'W(x,p') = |^(z)| 2 , (12) 



exactly as in the hydrodynamic equations, Eqs. (§,0). The result of Eq. (|Tl]) is outlined in Appendix [A]. 
If we take only the zeroth term in the infinite sum, we do indeed obtain a kind of Liouville equation 

— W = —H —W - —H — 

dt dx p dp dp 9 dx 

where 

„2 



W = — H P7T W ~ 7T H p— W ' ( 13 ) 



H p = ^ + V(x,t)+up, (14) 



V_ 

2m 

i.e. there is an additional "potential" proportional to the density of the distribution in position space. This can be 
interpreted as a large number of classical particles initially placed in phase space according to some kind of distribution 
function and interacting repulsively with one another, i.e. as a kind of non-ideal gas. If u is large we would generally 
expect large numbers of such particles concentrated heavily in some cell in position space to tend to drive one another 
apart, meaning that large values of p should in the long term be heavily disfavoured. 
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B. Hydrodynamics Related to Wigner Function Dynamics 



Hydrodynamic equations can also be derived from the equation of motion for the Wigner function Eq. (11), and 
if one expects the hydrodynamic equations to describe a semiclassical limit of the Gross-Pitaevskii equation, this 
should be consistent with the semiclassical limit described by the Liouville-like equation of Eq. (|l^) . In this section 
we conclusively show this not to be the case, and explain why this is so. 

In terms of the Wigner function, P is defined by: 

/oo 
dppW{x,p), (15) 
- OO 

where p(x) has already been defined by Eq. (|l2|). P is thus seen to be simply the first order momentum moment of 
the Wigner function. It turns out to be useful to define higher order moments as well: 

/oo 
dpp n W{x,p). (16) 
-oo 

The derivation of the equation of motion for p is carried out in Appendix [b|, and is exactly the continuity equation 
of Eq. (^) , correct to all orders in h. The equation of motion for P, again to all orders in h, turns out to be 

?rP = -i-W(x,t) + up] - — P 2 (x) + —^-( P P) 
at ox pm pm ox 
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where o~ 2 {x) — P2(x) — P(x) 2 is the variance of the Wigner function in p at a given point in x. 

Except for the term involving cr^, Eq. ( p"7| ) is identical to the hydrodynamic equation Eq. (|7]). However, it can 
be seen that Eqs. (§,|l^) do not form a closed system, as the equation of motion for P(x) refers to the higher order 
moment P2{x). There is in fact, as shown in Appendix [b|, an infinite chain of differential equations for the moments 
Pn{x) @: 

^-Pn(x) = J-[pP(z)] - — ^-[ P P n+ i(x)} - nP n ^(x)^-[V(x,t) + up] 

at p ox pm ox ox 

^ f (H/2) 2s r)l d 2s+1 1 
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In each equation the quantum corrections are described by the sum, but there is also an infinite chain of classical 
corrections; the second term of Eq. ( |l8| ) refers to the higher order moment P„+i(x). To get the second hydrodynamic 
equations, Eq. (^), in closed form from Eq. (|l7|), we must additionally make the zeroth order moment approximation, 

P n (x)=P(x) n . (19) 

In Appendix [b| this is treated in a little more detail. 

In order to reach the "hydrodynamic limit" , it is necessary to kill off all the quantum corrections, but there is in 
fact a much more drastic approximation than only taking the limit h — » 0, as a whole chain of classical corrections 
must be abandoned at the same time. The reason for the failure of the hydrodynamic equations as a semiclassical 
limit can be seen by examining our initial reasoning more closely. This was based partly on a correspondence between 
the hydrodynamic limit of the linear Schrodinger equation and the equivalent Hamilton- Jacobi equation, however this 
also implicitly assumes that the interpretation of the quantum wavefunction tends to a classical point. The Liouville 
dynamics given by Eqs. (10 13) describe the motion of classical distributions. As has already been mentioned, in 



the case of no nonlinearity(u = 0) one can connect the two classical cases by considering a distribution of the form 
W = 5(x — Xo)8(p — po), but when one is considering a case where the dynamics are influenced by the density in 
position space p, this is clearly meaningless. 

The correct semiclassical limit described in terms of moment equations is thus described by the following system: 

Ft P ^ = ~ p^h pP ^ - nPn^iVM^ (21) 
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where we must include every value of n. All of this is accounted for in Eq. (|T|). It seems clear that Eq. ( |l3| ) is a 
simpler way of describing the correct classical limit, and is almost certainly easier to integrate numerically. 

The purpose of comparing a nonlinear Schrodinger equation with its semiclassical limit is that it explicitly removes 
the wave-like or quantum behaviour, allowing us to see what there is that is specifically "quantum" about the dynamics 
of the nonlinear Schrodinger equation under consideration. 



IV. MODEL 



A. The Delta-Kicked Harmonic Oscillator 



To gain insight into the general problem, it is useful to take a simple test system, which is (a) accessible exper- 
imentally, and (b) amenable to numerical attack. The system chosen is the one dimensional delta-kicked harmonic 
oscillator, which has been studied both classically fl4|-|l7|| and quantum mechanically [[l8|-f2l| . The total potential for 
the classical Hamiltonian consists of a standard harmonic potential perturbed by a time dependent kicking potential: 

2 2 00 

V(x,t) = mUJ 2 X + Kcos(kx) S{t-nr), (22) 

n— — 00 

where x is the position, m is the particle mass, u the harmonic frequency, K the kick strength, k the wavenumber, 
and t the time interval between kicks. 



B. Scaling 



In our model, there are two basic parameters: the kick strength K, and the strength of the nonlinearity u. Ad 



ditionally there is h, which we have expanded around in Section [II. The parameters K and u need to be rescaled 
so that they remain equivalent in different regimes, as determined by a scaling parameter which takes the place of h. 
In the case of the delta-kicked harmonic oscillator there is a natural dimensionless scaling parameter, which is 77, the 
Lamb-Dicke parameter. 



" = fc fc (23) 

It should also be pointed out that r\ is a real physical magnitude, which really can be adjusted in the laboratory, 
unlike H. We call the dimensionless kicking strength k, and the dimensionless nonlinearity strength v. 

Kk 2 ^ 
\f2muj 1 ' 

uk 3 , , 

v = —= . 25 

2\/2mw 2 

It is shown in Appendix |c| that k and v have an equivalent effect on the overall dynamics for any value of r\. 

If, as is often the case when the trapping potential is harmonic, the Gross-Pitaevskii equation has been rescaled in 
terms of harmonic coordinates [xh = \J muj/hx, ph = p/V rnhuj), then it can be written in terms of these dimensionless 
parameters as: 



d 1 d 2 



u 



^^-^^^ + ^( a; ^V + ^MV (26) 
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V(x h ,t h ) = y + -j^S cos(V2nx) " nr '0- (27) 



n— — 00 



The wavefunctions have been rescaled so that they are properly normalized with respect to the harmonic position 
coordinate, and the time evolution is with respect to the dimensionless time th = tot. It is this form of the Gross- 
Pitaevskii equation that we use in our numerical simulations. 
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V. MODEL PHASE SPACE DYNAMICS 



A. Classical Point Dynamics 

The dynamics of a classical point particle in a delta-kicked harmonic potential have been described fairly extensively 
elsewhere fl^-|l7| . Briefly, we choose a value for th ■ For a given r there is only one free parameter which affects the 
phase space dynamics: k. There is a resonance condition = 2nr/q (rjq is a positive rational, where q > 2), whereby 
there are interconnecting channels of chaotic dynamics in the phase space [ fL6| , fl7| , the thickness of which depends on 
the kick strength k |[(|. For k not too large, these form an Arnol'd stochastic web which spreads through all of phase 
space, and has a characteristic q symmetry. For large k, one observes global chaos. Note that Arnol'd diffusion |22j ] 
can occur in systems of less than two dimensions when the conditions for the KAM (Kolmogorov, Arnol'd, Moser) 
theorem [p][2^] are not fulfilled, as is the case here JlB-Q. 

Here [and also in the following numerical work on the Gross-Pitaevskii equation Eq. (0) and Liouville equation 
Eq. (|l3|)] we consider the case where t/j = 2ir/6 and k = 1 . The scaled position and momentum are defined as 

x=^=T]x h , (28) 

k ' P =VPh- (29) 



\f2mi 

These scaled variables are chosen so that the phase space dynamics of a classical point particle described in terms 
of them are affected only by k and tv As can be seen, they correspond exactly to the scaled harmonic position and 
momentum when T) = 1. 

It can be seen in Fig. || that the phase space, in this case having a 6 symmetry, consists of a stochastic web of 
chaotic dynamics, where an initial condition can spread throughout phase space, enclosing cells of stable dynamics. 
An trajectory initially inside one of these stable cells will generally be held in a ring of six cells, equidistant from the 



centre, for all time (with the exception of the particle initially in the central cell, where it stays) 14 



B. Gross-Pitaevskii Equation 



In this section and in Sec. VC, we always work with the harmonically scaled position Xh and momentum p^ and 
with the dimensionless time th- For the sake of brevity we omit the h subscript, and thus write these variables simply 
as x, p, and t (or r). 

We integrate numerically, using a split operator method, the Gross-Pitaevskii equation as given in Eq. (|6|) consid- 
ering only the harmonic potential for periods of time of length t, punctuated by the exact mapping 

<p(x, t+) = e-^cos^zVvW^ r j ( 30 ) 

which accounts for the effect of the instantaneous delta kicks. This was carried out for various values of v and 77, 
where n = 1 and r = 2it/6 in every case. 

We have calculated the time averaged Wigner function, by which we mean the average of all the Wigner functions 
determined just before each delta kick, for 100 kick periods. The initial wavefunctions are displaced ground states. 
That is, the ground state of the Gross-Pitaevskii equation is determined numerically, for each value of v. We then 
locate the centre of the wavefunction at a point which is in a regular or chaotic region of the the classical single 
particle phase space. "Unstable" initial wave-packets are centred at x = y/2n jr\ (harmonic units), and "stable" initial 
wave-packets at x ~ 2^/2tt/i]. The initial wavefunctions are thus centred exactly either in the middle of a cell in phase 
space, or in an area dominated by web dynamics. These displaced states are the natural equivalent of coherent states 
for a cubic nonlinear Schrodinger equation. Just like coherent states, the density profile keeps its shape in a simple 
harmonic potential as it oscillates back and forth. This oscillating excitation is the so-called Kohn mode [pif . 

Firstly we show, in Fig. (||), the case of no nonlinearity, for the sake of reference. In this case the initial conditions 
are simply coherent states. Note that because it is possible for the Wigner function to have negative values, the 
colour representing zero is in general different in each pscudocolour plot. Thus, in each plot there is a "background" 
colour, which represents zero, with a superimposed pattern made up of darker and lighter shades. Notice that for 
77 = 1, the unstable initial condition [Fig. |(a)] appears to move through phase space following the stochastic web, 

iig.gb 



whereas the stable initial condition [Fi; 



(b)] simply circles around phase space, as would an initial coherent state 



in a simple harmonic potential. The wavefunction is clearly somewhat deformed (in the case of a harmonic potential 



G 



we would see perfect circles) but is otherwise well localized and well behaved. In the case of i] — 2, one might be 
forgiven for thinking that whether the initial condition is ostensibly stable or not is of negligible importance. The 
fact that rj is larger has the effect that the phase space is smaller compared to the size of the initial wavefunction (as 
plotted here, using harmonic units), and also quantum corrections play a bigger role (see Appendix |c|) , leading to 
the "tunneling" seen in Fig. ||(d), through classically forbidden areas of phase space. This tunneling can take place 
because the eigenstates of the Floquet operator F describing the period from just before one kick to just before the 
next, 

p _ ^-i(x 2 +p 2 )r /2 e ~iK cosfv^rji)/^') 2 {3^) 

are highly delocalized |lg|-|2"o[] , as is described in Appendix [d|. 

In Fig. |^ equivalent plots are shown when a nonlinearity of v = 0.1 is added to the Gross-Pitaevskii equation. It 
can be seen that this does not make very much difference to the phase space dynamics compared to no nonlinearity 
(Fig. |^), which is not really unexpected. 

When, as shown in Fig. |[ a nonlinearity of v = 1 is added to the Gross-Pitaevskii equation, it can be seen that 
this does make a difference. Intriguingly, given that the interaction potential is more strongly repulsive, the phase 
space dynamics appear to be more strongly localized. In the case of an unstable initial condition [Fig. |(a] and ||(c)] 
the web structure is noticeably reduced, and whereas in Fig. ^|(d) there was significant tunneling leading to a very 
delocalized phase space distribution, in Fig. ^(d) this has effectively disappeared. 

In Fig. [| where v = 10, this is even more marked. Where r\ = 1, in the case of an unstable initial condition 
[Fig. ^](a)], density seems to be concentrated around a "ring" in phase space, based around how far out in phase space 
the initial condition was. Where i] = 2 [Fig. 0(c,d)], whether the initial condition is ostensibly stable or unstable, we 
see only six symmetrically placed round blobs of density analogous to a coherent state in a harmonic potential. 



C. Liouville Equation 



Here we wish to investigate the semiclassical limit of the dynamics of the Gross-Pitaevskii equation with a delta- 
kicked harmonic oscillator potential [Eq. (p6|)]. The appropriate dynamics are described in general by Eq. (13). As 
with the Gross-Pitaevskii equation, in our case this can be carried out by considering only the harmonic potential for 
periods of time r, punctuated by an exact map describing the momentum kick. 

Equation (H|) can be qualitatively determined numerically by taking an ensemble of starting points from some 
desired distribution, using Hamilton's equations of motion to determine the trajectories, and using the numerically 
determined coarse-grained density for the overall potential governing the motion of the individual points. Obviously 
the coarse-grained density must be determined sufficiently frequently so that between times when it is determined, 
it does not change enough to have a very significant effect on the dynamics. This is in some sense analogous to the 
split-step method we have used to integrate the nonlinear Schrodinger equation, where as the time steps shrink to 
length zero, the approximate solution converges (in principle) to the exact solution. 

In each case the initial distribution is chosen by determining the ground state of the harmonic potential Gross- 
Pitaevskii equation (for appropriate v and 77), shifting it so that the centre of the wavefunction is at an unstable or 
stable fixed point (in the classical, single particle sense), calculating the Wigner function, and interpreting this as a 
classical probability distribution in x and p. The ground state Wigner function in the case of a harmonic potential is 
always strictly nonnegative, so one can always do this. 

Note that although 77 does not enter into the dynamics of Eq. (13) directly, by the above recipe it does enter by way 
of the choice of the initial condition, which affects the effective potential due to the distribution's density in position 
space, and so on. The time averaged density distribution plots in Figs. ^-|| are chosen to have initial conditions and 
scaling exactly equivalent to the time-averaged Wigner function plots shown in Figs. [l| [tj. 

In Fig. |?] we see the density distribution averaged over 100 kicks for the case where v = 0.1. The dynamics are 
essentially similar to those show in Fig. |^ for various single trajectories, and we observe a much lesser degree of 
distribution through phase space when compared to the full Gross-Pitaevskii equation (see Fig. |J). In particular 
we see no tunneling in Fig. [7(d), compared to Fig. ||(d). The dynamics in the cases of "unstable" initial conditions 
perhaps do not appear to be very strongly chaotic. Remember that only 100 kicks have been applied, and that in 
the case of the single particle classical delta-kicked harmonic oscillator, there are slow chaotic dynamics along the 
stochastic web jjTH16|, with an overall tendency to diffuse "outwards" in phase space. We have examined the case of 
100 kicks only in order to directly compare with the the numerically determined Gross-Pitaevskii dynamics. 

If we examine Fig. ||, which shows analogous dynamics to Fig. [?] for the case that v = 1, we observe some increased 
spreading out through phase space, still contained within the characteristic cells formed by the stochastic web in the 
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case of the stable initial condition for 77 = 1. In the case of rj = 2 The initial distribution seems too large for the cells, 
and even in the stable case there is some diffusion outwards through phase space. 

Finally we consider the case where v = 10, shown in Fig. ^. There is significant additional diffusion through phase 
space for the unstable initial condition, compared to the cases of v = 0.1 (Fig. |) and v = 1 (Fig. ||). Even for the 
supposedly stable initial condition there is some density which has found its way onto the stochastic web, and appears 
to be diffusing outward. Nevertheless, the basic structure of the single particle stochastic web appears to be retained. 

There thus appears to be a clear trend, where the larger the interaction parameter v, the greater the degree of 
diffusion outward through phase space, but nevertheless along routes typical for single particle dynamics. This has 
a simple explanation: when v is large and the distribution is highly localized, the distribution tends to push itself 
apart. After this initial explosion through phase space (actively encouraged in the unstable parts of phase space) the 
contribution by the density to the effective potential is small, and so the by now thinly spread distribution undergoes 
local dynamics equivalent to single noninteracting classical particles, chaotic or stable, depending on the location in 
phase space. 



D. Interpretation 



1. Overview 



The most interesting thing shown by these numerical experiments, is the conclusive demonstration that the local- 
ization observed in Figs. ^|, ra is due to interference effects, caused by terms of higher order in h in Eq. ( |i"l| ) (or more 
correctly, higher order in 77 , as shown in Appendix. |^). The intuitive picture of a stronger repulsive interaction 
driving the Wigner function/Liouville distribution apart, is fulfilled in the semiclassical limit, but breaks down when 
all "quantum" corrections are accounted for. 

The increasing degree of localization shown with increasing v in the Gross-Pitaevskii dynamics can also be qual- 
itatively explained. As is shown in Appendix in the case of linear Schrodinger equation dynamics, the Floquet 
eigenstates are highly delocalized, due to extra symmetries connected to the fact that the wavefunction is kicked 
exactly six times per oscillation period. The presence of delocalized eigenstates means that the wavefunction tends to 
spread throughout phase space with ease; along the stochastic web if the initial condition is in a classically unstable 
part of phase space, and possibly by tunneling from cell to cell (promoted by large 77) if the wavefunction is initially 
in a stable part of phase space. With increasing v this symmetry is more and more perturbed, to a point where this 
ability to spread freely through phase space is lost. Interference effects due to higher order terms of the density in 
Eq. (O), act to hold the wavefunction together, in contrast to the Liouville type dynamics described by Eq. (|l3|). 



2. Density in Position Space 



On this note it is instructive to look at the kinds of densities actually produced. We consider the final wavefunction, 
produced after 100 kicks, at a time just before a hypothetical 101st kick. In Fig. [To] we see plots of |</5(a;)| 2 for the 
case v = 0.1. Unsurprisingly for the unstable cases, and also for the stable case where 77 = 2, the states are highly 
delocalized in position space, with a great deal of fine structure. In Fig. [ll], this has substantially changed; the 
densities which were very complex are now much simplified, and even the stable initial condition for_n = 1 appears 
to have less structure when v = 1 compared to v — 0.1. When v is increased to 10, as shown in Fig. [l^, there is still 
some structure to the densities where 77 = 1, wheras in the case where 77 = 2 there appears now to be none. 

Obviously much more radical change is induced for the case of 77 = 2 when increasing v. Bearing in mind that rj 2 is 
our effective H, it is clear from Eq. (11) that higher order derivatives in the effective potential V(x, t) + up(x) will be 
more strongly emphasized (see also Appendix O). Between kicks, the non-Liouville corrections are due only to p(x), 



Considering the cases of Figs. |12|(a) and |12j(b) in particular, one might ask what there is about these densities which 
seemingly so totally dominates the dynamics. We consider the initial state, which is simply a shifted ground state. 
The ground state of the Gross-Pitaevskii equation lie somewhere between the cases of a Gaussian (no nonlinearity) 
and the Thomas- Fermi limit which is essentially an inverted parabola (large nonlinearity). With regard to the 
parameters we have chosen to use, the degree of "Thomas-Fermi-ness" is proportional to v/i] 3 . In the Thomas- Fermi 
limit, there are no higher order derivatives of p. A Gaussian however, has an infinite number of derivatives. For 
Figs. |l^(a,b), ?j/?7 3 = 1.25 only. The initial state density is thus more Gaussian than Paraboloid, and the large value 
of the effective h ensures that corrections due to the inevitable higher order derivatives are substantial. 
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Briefly: the application of a kick scrambles the phase of the position representation of a wavefunction; instanta- 
neously the density in position space is unaffected however. When looking at Eq. (|lTJ) we see that corrections due to 
higher order derivatives of p will be emphasized for larger effective H, in our case rj 2 . The effect of these corrections 
appears to be a strong tendency for the shape of the wavefunction to be preserved. 

In this work, we have not really explored the regime of very large nonlinearities. In view of the fact that in the 
Thomas- Fermi limit for the ground state there are no corrections to the Liouville-like equation of Eq. (|l3|), it is 
possible that the kind of very pronounced localization observed for the case of 77 = 2 might again be suppressed for 
much larger v. 



3. Density in Momentum Space 



For the sake of comparison, in Figs. ( |13| , [14 15) we show the corresponding momentum densities to the position 
densities of Figs, ([l^ jll],[l^) . The densities in position and momentum space essentially correspond, in that complex 
structure in one indicates complex structure in the other. This is not surprising, if we consider the kinds of Wigner 
functions displayed in Figs. fflEjJjl). 



VI. PHYSICAL MODEL: DRIVEN BOSE-EINSTEIN CONDENSATE 



A. Introduction 



A series of pioneering experiments investigating quantum chaos with atom-optical systems has been carried out by 
Raizen and co-workers mainly for a quantum realization of the delta-kicked rotor. We take a similar approach; 
a possible physical realization of the delta-kicked harmonic oscillator, consisting of a single trapped ion periodically 
driven by a laser, has been described in |26f| . This can in principle be readily extended to a periodically driven 
Bose-Einstein condensate. 



B. Single Particle 



We begin by regarding a single two level atom. In the x direction, it is trapped in a harmonic potential of frequency 
w, and driven time dependently by a laser field of Rabi frequency f2(t), wavenumber k, and frequency u>l- We disregard 
motional degrees of freedom in the y and z directions as being presently uninteresting, and arrive at the following 
Hamiltonian operator: 

H= h + + 2 {Wo(|e>(e| " l5><ffl) + cos (^)[^K Mit |e)<5l + H.c.]}. (32) 

In a rotating frame defined by 

U = exp[-WLt(\e)(e\-\g)(g\)/2], (33) 

and in the limit of large detuning |A| = \lol — ujq\ ^> |£2(i)|, |e) can be adiabatically eliminated to give, after 
transformation to an appropriate rotating frame, 

^=^ + ^ + ^[ooa(2**) + %>< 9 |. (34) 

The laser is periodically switched on and off, giving a series of short pulses, approximated by Gaussians: 

00 

n(t) 2 = n 2 e-('-" r » 2 , (35) 

n— — 00 

which approximate a series of delta kicks in the limit a — ► 0. Note also that we require a ^S> 1/A, otherwise the laser 
is too spectrally broad. Thus, we have finally 

H = f- + + 5£^[oos(2*f) + l]\g){ 9 \ ± 5(t - nr). (36) 







Because we assume that the atom is always in electronic state \g), the \g)(g\ operator can be effectively abandoned. 
The extra +1 simply adds a global phase, which can easily be accounted for, and so this can be further simplified to: 

A = £ + + cos(2*4) ± S(t - nr). (37) 

n=— oo 

This is exactly the Hamiltonian for the quantum delta-kicked Harmonic oscillator, except that we have cos(2fcx) 
instead of cos(fcr). As far as scaling is concerned, this means we must in turn consider rf — 2ry instead of r\ as the 
appropriate dimensionless parameter. 



C. Many Particles 

It is clear that if we consider a many particle system, then the above derivation is independent of any particle- 
particle interactions which do not change the internal states of the atoms. We thus consider the model Hamiltonian 
of a weakly interacting Bose gas, in second quantized form: 



H 



dx^(x) 



-— V 2 + V(x, t) + f tft(f)tf (£) 
2m 2 



(38) 



where i& is the particle field operator, g = 4:irfi 2 a s /m, and a s is the s-wave scattering length. We take V(x, t) to be 



V(x,t) = V(x,t)+'-^(y 2 +z% 
where the potential in the x direction is exactly that derived above, i.e. 



V(x,t) = 



8A 



cos(2fca;) 5(t — nr). 



(39) 



(40) 



We assume the radial frequency uj r to be very large compared to the axial frequency lu (cigar shaped trapping 
configuration), and thus assume that every particle is in the harmonic oscillator ground state in y and z. With this 
assumption we can integrate over y and z, reducing to a single dimension: 



H 



dx&(x) 



.^^ + V(x,t) + 9 -f^(x)i>(x) 
2m ox 2, 2 



*(a:), 



(41) 



where gid = mujg/2Trh = 2fkv r a s 



D. Asymptotic Expansion 

Using the particle number conserving formalism of Castin and Dum p^ j , we split the field operator ^ of the many 
particle system into a condensate part and a non-condensate part: 

$(x,t) =cp ex (x,t)a ipex (t)+64>(x,t), (42) 

where ip ex is the exact condensate wave function, and 5^ describes the non-condensate particles. Introducing the 
operator 

A ex (x,t) = -4=aJ, (t)8*(x,t), (43) 

Vn 

it is possible to make asymptotic expansions of A ex (x,t), <p ex (x,t), such that 

: A(2) + • • • , (44) 

? (2) + • • • , (45) 



A + 


J_a( x ) - 


i 




VN 


N 


ip + 




1 

1- — 


Vn 


N 
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where N is the total particle number operator. 

Thus, to lowest order, the condensate particles are described by <f(x). The time evolution of this can be shown to 
be given by the Gross-Pitaevskii equation |27j , which in our case is 

d ti 2 d 2 

ih—ip = + V(x, t)<p + Ng ld \ V \ 2 <p, (46) 

where N is the total number of particles. In turn, the non-condensate particles are described to lowest order by 
k(x,t). 

The Gross-Pitaevskii equation which we have arrived at in Eq. ( filf ) can be rewritten in terms of the dimensionless 
parameters rj', k, and v, as described in Sec. IV B, where 



V' = k\l— , (47) 



moj 
hk 2 a^]2VL 2 

2rawA 
8hNk 3 uj r a s 
\/2mw 2 



(48) 
(49) 



E. Non-Condensate Particles 

The mean number of the non-condensate particles is given by (6N) — (S^d^f), which to lowest order may be 
described by (A T A). In turn, At and A can be expanded as 

( A(a;, t) \ _ ? f u k (x, t)\ sTft ( v%{x, t) \ , . 

{ At( S| t))-L hk \ v k (x, t)) + ^ bk { «;(*, t) ) ■ (50) 

which gives rise to the following equation describing the mean number of non-condensate particles to lowest order in 
the perturbation expansion: 

oo 

(6N(t)) = J2(blh)(u k {t)\u k (t)) + (b{b k + l)(v k {t)\v k (t)). (51) 
fc=i 

The b k are time-independent ^Z7§ . We see that the time-dependence of Eq. ( |5l| ) is thus contained completely within 
(u k \u k ), (v k \v k ). A system initially prepared at temperature T has (blb k ) — [exp(E k / ksT)]^ 1 , and so, if we take the 
limit T — > 0, we get 

oo 

{SN{t)) =^{v k {t)\v k (t)). (52) 
fc=i 

We thus wish to study the dynamics of \v k (t)} to get some idea of the change in the number of non-condensate particles, 
in an analogous fashion to the work of Castin and Dum when investigating the behaviour of a condensate held in a 
time dependent isotropic harmonic potential p8|| . Note that because the Gross-Pitaevskii equation is nonlinear, it is 
possible to have chaos in the sense of exponential sensitivity to initial conditions within the Hilbert space. If this is 
the case, the above estimate of (SN(t)) will grow automatically, due to the fact that this estimate is essentially from a 



linearization around the Gross-Pitaevskii solution |28|. Thus the rate of growth of this estimate of (SN(t)) is similar 



to the Lyapunov exponent for the divergence of trajectories in phase space for discrete classical systems. 
The dynamics of the \u k (t)} and \v k (t)} are given by 

iu ± ( \u k (t)) \ _ rM ( \u k (t)) 

where 



in jt US = £ (HS (53) 



r(f) = / H G p(t)+Ng 1A Q(t)\(p(x,t)\ 2 Q(t) Ng ld Q(t)cp(x,t) 2 Q*(t) « . 

' ' 1 Ng ld Q*(tMx,t)* 2 Q(t) -H GP -Ng ld Q*(t)\if(x,t)\ 2 Q*(t) > 
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and where we have defined the Gross-Pitaevskii "Hamiltonian" , 

a 2 

H GP (t) = ^ + Vfr *) + ^(i,i)| 2 - £(f). (55) 

The phase factor £(t) is equal to the ground state chemical potential /i when (p(x,t) is the Gross-Pitaevskii equation 
ground state, for a harmonic potential. The projection operators Q, Q* are given by 

= 1-1^1, (56) 
Q* = 1-|*>WI> (57) 
where \ip*) is defined by {x\<p*) = <p*(x) = (f\x). 



F. Dynamics of (SN(t)) 

To determine how (SN(t)) changes over time we need to determine the dynamics of \vk(t)), which are coupled to 
the dynamics of |«fc(i)) through Eq. (|53|). We thus need to integrate Eq. (|53|), and to integrate Eq. (|53|), we need as 
initial conditions \uk(0)), \V).(0))- 

The initial conditions \v,k(0)), \i>k(0)}, for ip(x) in the ground state for a harmonic potential are determined by 
diagonalizing C where <p(x, t) is chosen to correspond to the Gross-Pitaevskii equation ground state, for a harmonic 
potential, and = fi. For this we need to determine the ground state condensate wavefunction ip(x) and the ground 
state chemical potential /i. This is achieved by propagating the Gross-Pitaevskii equation in imaginary time, where 
we use a split-operator method. 

We then determine C in the position representation where tp(x, t) is the previously determined ground state and 
£(t) = [i. We use a Fourier grid [ p9f to describe p 2 in the position representation. We then diagonalize C numerically, 
and gain as the resultant set of eigenvectors 

with eigenvalues {-Efej — -Efe, 0, 0}, respectively ^7j. These eigenvectors must be properly normalized ]27[j , so that 

dxu* k (x)uk' (x) — / dxvl(x)vk' (x) = Skk' (59) 



Our initial condition for the Gross-Pitaevskii equation is in general a shifted ground state, that is, we take the ground 
state wavefunction, and instantaneously translate it in position space, otherwise altering nothing. Physically, this 
could be achieved by almost instantaneously translating the centre of the harmonic potential, so that x 2 — > (x — a) 2 . 
Instantaneously, this would leave the Gross-Pitaevskii wavefunction and the Uk(x), Vk(x) modes unchanged. If we 
then re-express everything in terms of x 1 = x — a, we end up with the same equations in terms of x' as we had initially 
in terms of x, but the wave functions are transformed: {(p(x),Uk(x), Vk(x)} — > {fix' + a), Uk(x' + a), Vk(x' + a)}. 

Thus, if the initial Gross-Pitaevskii wavefunction is simply a shifted ground state, then the appropriate initial Uk, Vk 
are correspondingly shifted from those determined from C for the ground state condensate wavefunction. This set of 
initial conditions is in fact somewhat special; as previously mentioned, the density profile of ip(x) remains unchanged 
as it oscillates back and forth (without kicks), the same is also true of Uk(x) and Vk(x). 

Once we have the initial conditions we can start integrating Eq. (|53|). 



G. Numerical Results 



We integrated numerically Eq. ( p3[) for the first fifteen Uk(x), Vk(x) pairs over the time span of 100 kicks, using 
a split operator method described in some detail in Appendix parallel to numerical integration of the Gross- 
Pitaevskii equation, also using a split operator method. Just before each kick each of the inner products (vk\vk) were 
determined, which are plotted against time in Figs. [l€-19, for various parameter regimes we have already investigated 



the Gross-Pitaevskii dynamics of. The "stable" and "unstable" initial conditions referred to are those of the initial 
Gross-Pitaevskii wavefunction [which in turn determines the initial conditions of each of the Uk (x), Vk(x) modes], 



and are exactly those taken in the integrations of the Gross-Pitaevskii equation described in Sec. II A. To reiterate 
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the data presented in the plots in this section correspond exactly to the phase space plots presented in Sec. Jlj_A| for 
the appropriate values of v and rf , with regards to the initial condition. Figs. |l6| , |l7| correspond to Figs. and 
Figs. |lq JTa co rrespond to Figs. |||. 

In Fig. (|16D, where v[ — 1 and v = 1, we see a marked difference between the "stable" and "unstable" cases. In 
unstable case we see much greater growth of the (vk\vk)- Interestingly, the k — 1 mode in the stable case does not on 
average seem to grow at all, instead undergoing quasiregular oscillations in time. The leading terms arc also different; 
k = 1 for the unstable case, and k = 2 in the unstable case. 

Compared to Fig. |l6|, the "stable" and "unstable" cases shown in Fig. [n] (where the only difference is that r)' = 2), 
appear comparatively similar. In particular there does not seem to be a great deal more growth of the (vk\vk) in the 
unstable case when compared to the stable case. 

We see the same pattern repeated in Figs, [l^ and |l^, where v is now 10. In Fig. [l8| the (vk\vk) very rapidly grow 
in the unstable case when compared to the stable case, whereas in Fig. |l9|, where 77' = 2, the difference is not nearly 
so marked (and in any case the grow th of the (vk\vk) is generally less). This reflects in some sense the observed 
Wigner function dynamics in Sec. II A , where there does not seem to be such a strong qualitative difference between 
the "unstable" and "stable" cases where 77' = 2 for any value of v, in contrast to the cases where 77' = 1. One 
should bear in mind that although the dimensionless nonlinearity strength v is the same in both Figs. |l| and [l?], 
the actual repulsive interaction Nu\d is proportional to v/rj' 3 . One might argue then that one would expect that 
there is generally less depletion from the wavefunction described by the Gross-Pitaevskii equation. The evolution 
of ip(x,t) is also important however: v/rj 13 — 1 where v — 1 and r\' = 1 is not that different from v/r)' 3 = 1.25 
where v = 10 and 7/ = 2, but the evolutions of the (vk\vk) are. There appears to be some correspondence between 
the Gross-Pitaevskii phase space dynamics shown in Figs. and the evolutions of the (vk\vk), in that when there 
is a significant difference between the "stable" and "unstable" cases, this shows up in the dynamics of the (vk\vk) 
corresponding to these different cases. Also a more "smooth" phase space plot (as for 77' — 2 compared to 77' = 1 in 
Figs. appears to correspond to a more "smooth" evolution of the (vk\vk) (Figs. p"7 ,|l9| compared with Figs. |l6| , p"8| ). 
As the equation describing the time evolution of the \uk),\vk) pairs is essentially the same as that describing the 
evolution of linearized orthogonal perturbations of the Gross-Pitaevskii wavefunction ||27| , this is not unexpected. 



H. Comparison with Experimental Parameters 

We first examine our best estimate for (SN(t)}, which is J2k~i( v k(t)\ v k(t)} , where t is expressed as the number of 
kicks. In Fig. ^ this is plotted for each case where v — 1 against the number of kicks, and in Fig. |l| for v = 10. 
Interestingly, for v = 1 and 77' = 2, total growth appears to be almost exactly linear in time, after a short buildup 
period; as noted before, growth does not appear to be that different when comparing the "stable" and "unstable" 
cases. For 77' = 1 however, there is a clear and substantial difference between the two cases. 

When v is increased to 10, as shown in Fig. |l], growth becomes more erratic. We see that for the "unstable" case 
where rf = 1, X)fcLi( w fc|' l; fc) ends up being very large, making it unlikely that an experiment for this parameter regime 
would follow Gross-Pitaevskii dynamics. The general pattern observed in Fig. |2^ is repeated here, but with larger 
numbers. Note however, that the beginnings of a clear differentiation between the degree of growth for the "stable" 
and "unstable" cases when 77' = 2 appear to be occurring; in both cases growth is certainly not linear with time. 

Overall, our results can be interpreted as similar to those obtained in p8| for the case of a time dependent harmonic 
potential. When one would expect classical chaotic behaviour, one observes rapid growth of the (vk\vk)- 

To examine the behaviour of a possible experimental realization of this scheme, we consider Rubidium 87, which 
has an s wave scattering length of a s = 5.1 x 10~ 9 m Q, and Sodium 23 (a s — 2.75 x 10~ 9 to) [[n). Substituting 
Eq. ( fl7j ) into Eq. (^), we can rewrite v, so that 

v = J ^-2Nuj r a s r]' 3 (60) 
V nw 

is expressed in terms of 77', which is more convenient for our purposes. Using Eq. (|60|), we get as a general relation 
for the number of particles N = X^/uj/uj ri where 

A= vl^ (61) 

The values of A in units of s -1 / 2 for the parameter regimes we have investigated are summarized in Table [j]. 

We let oj r = 10w, remembering that we should have Lu r significantly bigger than u>, we take this to be a reasonable 
minimum, bearing in mind that the values of the harmonic potential ground state chemical potential \i lie between 
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0.55 and 3.11 in units of tko, as shown in Table |. We then get N — is/^/uJ^, where v = A^/l/10. Numerical values 
for v in units of s" 1 / 2 , where u> r = 10u> are also displayed in Table. B. In principle this leaves us one free parameter 
to tweak; the smaller the radial frequency, the larger N can be, and the less significant the effect of the growth of 
the number of particles not described by the Gross-Pitaevskii equation. This would mean that we could reasonably 
expect to describe the dynamics of the particles largely with the Gross-Pitaevskii equation, with small corrections 
accounted for by Eq. (p3|). 

In practice trapping frequencies for alkali atoms such as Rubidium and Sodium lie between about 1 and 100 
Hertz. The growth of J2kLi( v k\ v k) in the "unstable" case where v = 10, rf = 1 is thus far too high for this simplest 
interpretation of the real dynamics. The cases where if — 2 look more promising, and here in fact the interesting effect 
of nonlinearity induced localization within phase space of the Gross-Pitaevskii wavefunction is even more pronounced. 
Also note that even for a small nonlinearity of v = 1, there is still a pronounced difference in the Gross-Pitaevskii 
equation phase space dynamics (see Fig. ||) compared to the case where there is no nonlinearity (Fig. for both 
if = 1 and if — 2, and here the numbers also seem more promising for the nonlinearity induced localizing effect to be 
observed, corresponding to our numerical integrations of the Gross-Pitaevskii equation. 



VII. CONCLUSIONS 



We have derived explicitly an appropriate semiclassical limit for a general cubic nonlinear Schrodinger equation, or 
Gross-Pitaevskii equation, and find it to be a Liouville type equation, with a term involving the density in position 
space. We have shown how and why this differs from the hydrodynamic limit of the Gross-Pitaevskii equation. 
In particular, this derivation shows how an eccentric wavefunction ip(x) can produce large deviations from this 
semiclassical limit, through higher order corrections involving derivatives of the density p{x) = \(p(x)\ 2 , in addition 
to effects due to an unusual potential. We have investigated numerically a simple test system, the one-dimensional 
delta-kicked harmonic oscillator, studying the dynamics of the Gross-Pitaevskii equation and the appropriate Liouville 
type equation. We have found for moderate nonlinearity strengths that there is a localization effect explicitly due to 
interferences caused by the nonlinearity. We have outlined a possible experimental implementation of such a system 
in a Bose-Einstein condensate experiment, and have investigated numerically to what degree the Gross-Pitaevskii 
equation describes correctly the dynamics of the bulk of the particles for certain test cases. From this we have 
determined a lowest order estimate for the growth in the number of non-condensate particles. We have found that for 
this system this depends strongly on the parameter regime of rf and v under study, and that this seems to correspond 
to the kinds of phase space dynamics observed in the Gross-Pitaevskii equation. We have compared the numbers 
obtained with realistic experimental parameters for condensates formed from sodium or rubidium atoms. 



ACKNOWLEDGEMENTS 



We thank J. R. Anglin, for helping clear up a number of points on the work in Sec. IIIB , M. G. Raizen, Th. Busch, 
and K. M. Gheri, for discussions, and D. A. Steck for bringing reference [[l9| to our attention. We also thank the 
Austrian Science Foundation, and the European Union TMR network ERBFMRX-CT96-0002. 



APPENDIX A: DERIVATION OF WIGNER FUNCTION DYNAMICS 



1. Definitions 



Defining the Wigner function for a pure state as 



1 f°° 

W(x,p) = — / dTe-^/Vte - r/2)ip(x + r/2), (Al) 

27Tft /_„ 



we take the time derivative 



o o o 

g-W(x,p) = -W(x,p) SP + -W(x,p) NL , (A2) 

where we have split up the differential equation into a part which is governed by the single particle linear dynamics 
(SP), and a part which is governed by the nonlinearity (NL). 
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2. Single- particle dynamics 



The single particle dynamics are described by: 



dTe -irp/h 



(<p\H\x-T/2)(x + T/2\<p) - (<p\x-t/2)(x + t/2\H\<p) 



(A3) 



The expansion we desire is exactly that used by Zurek and Paz in investigating the quantum-classical boundary 0, 
and is based on work originally carried out by Moyal [B2| and Wigner pq] : 



§-W(x, P ) SP =± ( U 



h\ 2s d 2s+1 d 2s+1 a a 

H——W-—H—W. 



^(2s + l)!V2/ dx 2s+1 dp 2s+1 dp dx 



(A4) 



3. Nonlinear Dynamics 

For a simple cubic nonlincarity u|(^| 2 (/9, we can express dW(x,p)^/dt as 

dp' [W(x - t/2, P ') - W(x + r/2,p')] / dp"e lTp " /h W{x,p") I 



d iu 



We expand W(x — r/2,p') — W(x + r/2,p') as a McLaurin series: 



d in ^ (1/2)^+1 r- 



tt/i 2 ^ (2s + 1)! 



(9 



2s+l 



Using the chain rule and Fourier's integral theorem, we arrive at 
d , . u ^ (-ft/2i) 2s+1 9 2s+1 



tt/i 2 ^ (2s + 1)! dx 2 ^ 1 

s— 



dp'W(x,p') 



(A5) 



rfp' fe2s+1 VF(x,p') / dre- lTp/ri / dp"T 2s+1 e lrp " /h W(x,p"). (A6) 



(A7) 



4. Combined Result 



Combining Eqs. (A4) and (A7), we get the Wigner function dynamics to all orders in h of the cubic nonlinear 
Schrodingcr equation with arbitrary potential, in one dimension 







dt ^ 



(-l) s 



(2s + 1)! \2 J dx 



fc\2s g2s+1 Q2s+1 g Q 



Qp2s+1 Qp Q x 



which has as its semiclassical limit (ft — > 0) a Liouvillc-likc equation: 

d d d d d 

° W = 1L[H + up ] ° W _ ° H £-W, 
at ox op op ox 



(A8) 



(A9) 



where p is the Wigner function integrated over p, as defined in Eq. (|l2|). This derivation can be easily generalized for 
other nonlincaritics and to two and three dimensions. 



APPENDIX B: RE-DERIVATION OF THE HYDRODYNAMIC EQUATIONS 

1. Definitions 



The density p has already been defined in terms of the Wigner function by Eq. (|12|). The quantity P is defined in 
terms of the Wigner function as 



pP 



dppW. 



(Bl) 
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2. Regaining the First Hydrodynamic Equation 



The equation of motion for p is given by 



2s Q2s+1 



dt P ~ ^ (2s + 1)! \ 2) dx 2s + 

s—0 x / \ / 



T [H + up] / dp 



dp 



2s+l 



w 



dp—W—H. 

ox op 



(B2) 



Due to the fact that W(x,p) and all of its derivatives are equal to zero at x = ±oo, something we make frequent use 
of, this simplifies to the continuity equation 



-P = -- — { P P) 
dt m dx ' 



(B3) 



using the definition of Eq. (Bl) 



3. Equations for Higher Order Moments 



We now turn to the equation of motion for P. We have, from Eq. (Bl) 



q i r°° ( 00 C— i) s / h\ 2s d 2s+1 q2s+i q \ PQ 

—P = - dpply-^ — - n 9 -ii [V(x,t)+up] n „ ^ W-——W } + — TripP)- 
dt pj-oo 1 ^ (2s + l)! V 2 / 9x 2s+lL v ' ; Hl dp 2s+> ">>>,■ f :„,,,), ' 



to 9x 



/9TO 9x 



(B4) 



The integral of the Wigner function over p, dppd 2s+1 W/dp 2s+1 , is equal to p when s = 0, and is otherwise equal 
to zero. We therefore have 



d_ 

dt 



9 rTr/ v , 19 

Q- x [V( X ,t) +UP }--^ 



d PP 2 w) + —^-{ P P). 
pm ox 



(B5) 



Clearly Eq. (|B3| ) and Eq. ( |B4| ) do not form a closed system of equations, due to the presence of the second order 
moment P^ix), where 



Pn(x) 



p(x) 



dpp n W(x,p). 



It is relatively simple to derive a chain of equations of motion for all P n (x): 

P or 



d 1 f°° d 

-P n {x) = - / d PP n -W 
dt p J_ OQ dt 



Substituting in Eqs. ( |A8| , [B3| ), we get as the general form: 

1 d 



d Pn ( x) = ^MJL [pP{x)] _ J_|_[p P (:b )] - n P n ^(x)^-[V(x,t) +up] 
ot pm ox pm ox dx 

- E { ( 2s + i) ! V-( 3 + i) ]! p »-(^)W^Tr[^ < ) + «d} • 



(B6) 



(B7) 



(B8) 



The system of equations Eqs . (|B3| , [b"8| ) , where n ranges from 1 to oo, thus describes the full dynamics of the Gross- 
Pitaevskii equation, Eq. (Q) |13| . 



4. Regaining the Second Hydrodynamic Equation 



We consider a set of solutions of the moments where P n {x) = P(x) n . Taking Eq. ( |Bq ) and setting h = 0, i.e. 
ignoring all quantum corrections, we substitute this solution in, which after differentiation results in: 



nP{xY-^P{x) 



.^Ip^-nPixT-^V^+upi 



(B9) 
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where we can immediately carry out cancellations, to finally arrive at 



d_ 

dt 



P(x) 



d_ 

dx 



P{xf 
2m 



V(x, t) + up 



(BIO) 



which is the second hydrodynamic equation, Eq. (Q). Thus hydrodynamic equations describing dynamics in the 
hydrodynamic limit [§||l0|] are valid whenever h — > and P n (x) — P(x) n . This condition can be expressed in terms of 
Liouvillc distributions as 



dpp n (x)W(x,p) 



dppW(x,p) 



(Bll) 



which is in general fulfilled for W(x,p) = p(x)5[p — po(x)], where po{x) is some single valued function of x. 



APPENDIX C: MORE SCALING 



As dimensionless parameters we have m k, and v, defined in Eqs. ([25 ]2~4 ,^5| ) , respectively. We have as dimensionless 
coordinate and canonically conjugate momentum the variables of Eqs. (2^2^), and use the dimensionless time th = tot. 
Using this, we can write the dimensionless single particle Hamiltonian function as 



~2 
P 



H = ij + V(x,t h ) 



_ 2 txJ 

V(x, t h ) = — + — cos(\/2£) ^2 - nl ~h) 

* n— — 00 



the Gross-Pitaevskii equation, Eq. ([!]), as 

. d _ 



2 dx 2 

and the equation of motion for the Wigner function as 



Tj d^ 1 1 ~ ~ V 



dt h ^ 



H + vp 



d 2 ^ 1 . d ~ d ~ 
_ w - —H—W. 

Qp2s+1 Qp Q% 



{2s + 1)! V 2 / dx 2 ^ 1 

The wavefunction, Wigner function, and density, have been rescaled so that they are properly normalized: 

(p = \J V2/ktp, 



k 2 



p 



dpW. 



(CI) 
(C2) 

(C3) 
(C4) 

(C5) 
(C6) 
(C7) 



In the expansion shown in Eq. (C4) it can clearly be seen that if r\ is varied, then this is completely independent of all 
other rescaled quantities. We see thus that rj 2 is an appropriate expansion parameter, and that the other dimensionless 
parameters k and v are correctly scaled to be independent of the expansion parameter. If one takes only the zero 
order term in the sum, rj drops out completely. 



APPENDIX D: CRYSTAL SYMMETRY AND NON-LOCALIZATION 

1. Classical Background 

Consider the classical delta kicked harmonic oscillator described in Eq. |Cl]. The symmetry properties of this system 
have been extensively investigated by Zaslavsky and coworkers |l4|-|l^]; we recapitulate some of this to provide context. 
One can determine a kick to kick mapping terms of a = (x + ip)/^/2: 
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Oi n +l 



a n + i—= sin(a„ + a*) 
v2 



If ujt = 2nr/q, then we can write the mapping after q kicks as 

9-1 



—= sm{a n+k + a* n+k )e 
^ l fc=0 



%2-nkr J q 



(Dl) 



(D2) 



Keeping terms in k up to first order only, we observe an approximate rotational q symmetry in phase space ]l5| , pT[ ; 
if we substitute a n with /3 n = a n e l27Tl / q , ( £ Z, we end up with f3 n + q — a n + q e l27Tl ^ q . There can also be a translational 
symmetry in phase space, i.e. (3 n — a n + 7 [3 n +q = otn+q + 7,7 € C. Note that it is only possible to combine a 
rotational q symmetry with translational symmetry when q (z q c = {1, 2, 3, 4, 6} 
Translational symmetry demands 



9-1 

sin(a 

k=0 



n+j r a n+j. 



a i2-Kkr I q 



9-1 



(D3) 



k=0 



which in turn implies I3 n +j + /3* +fe = a n +k + a n+k + 27r/fe; V fc, Zfc € Z. Thus, Eq. (D2) for f3 n +q can be simplified to 



9-1 



/?n+<2 = a n + 70 + i-= sin(a„ +fc + a* +fe + 7 fe + 7 / *)e l27rfcr/9 , 

fc=0 



(D4) 



where 7^ = je j27rfcr /9 . The condition for translational symmetry is thus reduced to 7/. + 7^ = 2ttIk, which implies 

(D5) 



7 — 7 

ifc = Iq cos(2nkr/q) — i sin(27rfcr/g). 

27T 



If we now let k± — q/2 ± m or (g ± m)/2, depending on whether or not (7 is even, we get 



cos(2 7 rfc + r/( Z ) = lk ++ lk ~ e Q, 

7 — 7* 

i sin(27rfc + r/g) = — Zfc, G Z. 



(D6) 
(D7) 



This implies that cos(27r/g) G Q, and it is known that this can only be true if q € g c = {1, 2, 3, 4, 6} J35|. This directly 
implies that cos(27rA:r / q) € Q V fc,r € Z. There is is thus an e:raci translational or crystal symmetry in ph ase spa ce, 
for q G q c only. There are an infinite number of values of 7 for which this applies, determinable from Eqs. (D6 D7). 



2. Quantum Expression 

Broadly following the treatment of Borgonovi and Rebuzzini |l9|] , we consider the unitary displacement operator 

D{a) = e aa ~ a a — e % ^" x ~^ pfj]. The operators a 1 ' and a are the quantum harmonic oscillator creation and 
annihilation operators, and the operators x, p, are scaled in harmonic units. The displacement operator acting on 
a wavefunction is a quantum analogue to translating a classical point particle in phase space. We now consider the 
Floquet operator F — e - I ( !lt <'+i/ 2 ) w '' e - !K ™ s [i(''+« t )]/v / 2') 2 anc j determine the commutation properties of it with the 
displacement operator. 

Using elementary properties of coherent states (3(| , it can be seen that 

9-1 

D(a)F q = JJ | e -i(a t a+l/2)27rr/ 9e -iKcos[i ) (a+a-a fc -a^)]/V2r ) 2 I ^(a), (D8) 
k=0 



where at = ae l2 * kr l q . The product of Floquet operators F q corresponds to the mapping of Eq. (D2) which we used 
to investigate classical symmetry properties. 

Thus, D(a) commu tes with F q if r/(ak + cv* k ) = V^rj^k = 2nlk, h G Z V k. Using this we arrive at, similarly to 



the derivation of Eq. (D5) 
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l k = l Q cos(2Trkr/q) - " '— sm{2nkr/q). (D9) 

V27T 

Analogously to the classical case, D{a) commutes with F q if and only if q G q c . This implies that for q S q c , the 
eigenstates of F q are invariant under certain displacements, of which there are an infinite number, and are thus 
extended. Localization is not expected to take place, similarly to the case of quantum resonances in a delta-kicked 
rotor MM. 



APPENDIX E: INTEGRATION OF THE £ EQUATION. 

From [^Tj, we know that 

d ( !»*(!)) \_r! Mf)> 



K = £ wo" '■ (E1) 



and that the corresponding time evolution operator 



U{t) 



{ Q*(t)) UMt) { Q*(0)J- (E2) 



The operator Ugp (t) is the time evolution operator corresponding to Cqp (t) , given by 

r (A- ( V(x,t) + 2u\v(x,t)\ 2 +f/2m u<p(x,t) 2 \ . . 

L GHV-y ~u^{x,ty 2 -V(x,t)-2u\cp(x 1 t)\ 2 -p 2 /2m J ■ 

In our case, the potential is that of the delta-kicked harmonic oscillator. Integrating between kicks, we consider 
V(x) time independent. Note however, that £gp(^) is still in principle time dependent through ip(x,t). Thus, taking 
very small time steps At, the evolution is given approximately by 

\U k (t + At)) \ _ iC GP (t)At/H ( Wkit)) \ 

\V k (t + At)) )~ e {\V k (t))J- {M) 

The time evolution operator e" t£ GpW At / R can be split into position dependent and momentum dependent parts, and 
the time evolution was then determined using a split operator method, of which there are many variations f38| . We 
set \Uk{0)) = \uk(0)) and | Vfe (0)) = \vk(0)), and determined |itfe(t)) and \vt(t)) from \Uk{t)) and \Vk(t)) by projection, 
just before each kick . 

The effect of a kick is given by: 



(E5) 



In Sec. VI G, the procedure outlined above was used, in conjunction with numerical integration of the Gross-Pitaevskii 



equation, also by a split operator method with matching time steps. 
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FIG. 1. Schematic diagram of how nonlinear Schrodinger equations relate to other forms of dynamics under various limits. 
The parameter u represents the strength of the nonlinearity. 




FIG. 2. Poincare sections of the phase space dynamics of the classical delta-kicked harmonic oscillator, (a) Single unstable 
initial condition forming a stochastic web spreading through phase space, (b) Close up of the phase space, showing the closed 
curves characteristic of regular dynamics. In both cases 77, = 2tt/6, k = 1. 
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FIG. 3. Pseudocolour plot of time averaged Wigner functions when v = 0, i.e. linear Schrodinger equation dynamics, in the 
two cases of: r) = 1, for (a) unstable initial condition, (b) stable initial condition; r\ — 2, for (c) unstable initial condition, (d) 
stable initial condition. Position and momentum are scaled in harmonic units, and black means large and positive. 
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FIG. 4. As for Fig. H, where v = 0.1. 
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FIG. 5. As for Fig.H, where v — 1 
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FIG. 6. As for Fig. |[ where v = 10. 
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FIG. 7. Pseudocolour plot of time averaged distributions undergoing Liouville dynamics when v = 0.1, in the two cases of: 
r) = 1, for (a) unstable initial condition, (b) stable initial condition; r) = 2, for (a) unstable initial condition, (b) stable initial 
condition. Black means large and positive. 
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FIG. 9. As for Fig. | when v = 10. 




FIG. 10. Plots of |y(a;)| 2 after the application of 100 kicks and where v = 0.1, in the cases of: r\ = 1, for (a) unstable initial 
condition, (b) stable initial condition; and rj — 2, for (c) unstable initial condition, (d) stable initial condition. 
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FIG. 11. As for Fig. 0, but for v 




(b) 



0.2 



0' 

-20 -10 

(d) 



0.5 
0.4 
0.3 
0.2 
0.1 



-20 -10 



FIG. 12. As for Fig. [u| but for 
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FIG. 13. Plots ol |y(p)| 2 after the application of 100 kicks and where v = 0.1, in the cases of: 77 = 1, for (a) unstable initial 
condition, (b) stable initial condition; and rj = 2, for (c) unstable initial condition, (d) stable initial condition. 
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FIG. 14. As for Fig. jT|, but for v = 1. 
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FIG. 15. As for Fig. but for v = 10. 




FIG. 16. Semilog plot of change in (vk\vk) with respect to the number of kicks n, for k = 1,. . . , 15: fc — 1 solid line, fc = 2 
dotted line, fc = 3 dashed-dotted line, fc = 4 dashed line, fc = 5 circles, = 6 crosses, k — 7 pluses, k = 8 squares, = 9 
diamonds, fc = 10 downward pointing triangles, fc = 11 upward pointing triangles, fc = 12 left pointing triangles, fc = 13 right 
pointing triangles, fc = 14 pentagrams, fc = 15 hexagrams, where jj' = 1, and w = 1. (a) shows data for the "unstable" initial 
condition, where the leading term after 100 kicks is for k = 1, (b) shows data for the "stable" initial condition, where the 
leading term corresponds to fc = 2. 
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FIG. 19. As for Figjl|, except that Tj = 2, v = 10. In (a) the leading term is for k = 4, in (b) for k — 6. 
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FIG. 20. Plots of Ylk=i ( Vk \ Vk ) a g ams t the number of kicks n, where v = 1, in the two cases of: rj = 1, for (a) unstable 
initial condition, (b) stable initial condition; rj' = 2, (c) unstable initial condition, (d) stable initial condition. 
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Na 23 


Rb 87 


V 


i 


H 


A 


V 


A 


V 


1 


1 


0.87/kj 


9.55 x 10 3 s~ 1/2 


3.02 x 10 3 s~ 1/2 


2.65 x 10 3 .<T 1/2 


8.38 x 10 2 s~ 1/2 


2 


0.55hu 


1.19 x 10 3 s _i/2 


3.77 x 10 2 s _1/2 


3.31 x 10 2 s" i/2 


1.05 x 10 2 s" i/2 


10 


1 


3.Uhw 


9.55 x 10 4 s" 1/2 


3.02 x 10 4 s" 1/2 


2.65 x 10 4 s _1/2 


8.38 x 10 3 s" 1/2 


2 


0.95hw 


1.19 x 10 4 s _i/2 


3.77 x 10 3 s~ 1/2 


3.31 x 10 3 s" i/2 


1.05 x 10 3 s _i/2 



TABLE I. Values of A and v for Sodium 23 and Rubidium 87, when in the parameter regimes of v and r\ under study. Also 
displayed are the values of the numerically determined ground state chemical potential [i for the appropriate values of v and 
rf , in units of tko. 
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